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FINITE ELEMENT ANALYSIS OF NONLINEAR FORCED VIBRATIONS 
OF BEAMS AND RECTANGULAR PLATES 

Chuh Mei* and Kamolphan Decha-Umphait 
Old Dominion University, Norfolk, Virginia 

SUMMARY 

Slender beams and thin plate structures subjected to periodic loading 
are likely to encounter oscillations large in comparison to beam depth or 
plate thickness. The responses predicted using the small deflection linear 
structure theory are, therefore, no longer applicable. Nonlinear theory, 
taking into account the effects of large amplitude, must be employed. 

In this report, a finite element method is presented for the large 
amplitude vibrations of complex structures which can be modelled with beam 
and rectangular plate elements subjected to harmonic excitation. Both in- 
plane deformation and inertia are considered in the formulation. The method 
described would give more accurate results than the classic conti nuim 
approaches which neglect the effects of inplane inertia. Derivation of the 
harmonic force and nonlinear stiffness matrices for a beam and a rectangular 
plate element are presented. Solution procedures and convergence character- 
istics of the finite element method are described. Nonlinear response to 
uniform and concentrated harmonic loadings, and improved nonlinear free 
vibration results are presented for beams and rectangular plates of various 
boundary conditions. 


*Associate Professor, Mechanical Engineering and Mechanics 
tGraduate Student 


INTRODUCTION 


Linear models have had, and will continue to have, great importance in 
the analysis of real structural systems. There are, however, some dynamic 
systems which are inherently nonlinear and their behavior cannot be de- 
scribed fully through linear models alone. 

NASA is considering placing in orbit various large, flexible, and low 
density space structures. These structures, or some components of these 
structures, may be expected to have oscillations large in comparison to some 
corresponding dimension (beam depth or plate thickness). In this situation, 
the responses predicted using the small deflection linear structure theory 
may no longer be valid, therefore, nonlinear structure theory taking account 
the effects of large deflection may have to be employed. 

The first papers that dealt with nonlinear vibrations of beams and 
plates with immovable edge supports were those of Woinowsky-Kr ieger (ref. 1) 
and Herrmann (refs. 2 and 3). Herrmann presented a set of equations of 
motions that corresponds to the dynamic analogue of the von Karman plate 
theory. The nonlinear beam vibration was treated as a special case of plate 
vibration in reference 3. Although the basic equations had been thus 
established, the general solutions of the coupled, nonlinear, governing 
differential equations of motion are still not available owing to the 
mathematical difficulties involved. The classical continuum approach was to 
employ various approximate methods to beams and plates of simple geometrical 
shape and of simple restraint conditions along the edges. 

Nonlinear forced vibration of beams with different boundary conditions 
has been investigated by applying the Galerkin method (e.g. refs. 4 to 7) 
and the resulting equations were solved by the harmonic balance method. 
Other approaches to the problem, namely the finite element method combined 


with method of averaging, the form-function approximation method and the 
incremental time-space finite strip method were proposed by Busby and 
Weingarten (ref. 8), Lou and Sikarskie (ref. 9) and Cheung and Lau (ref. 
10), respectively. Further, a perturbation procedure based on the method of 
multiple time scales has been presented by Nayfeh et al. (refs. 11 and 12). 
Eisley (ref. 13) and Sathyamoorthy (ref. 14) have presented their 

comprehensive and excellent reviews on nonlinear analyses of beams 
concerning the classical methods. 

Nonlinear forced vibrations of circular and rectangular plates with 
various boundary conditions have also been investigated by using the 
Galerkin or Ritz method (refs. 4, 5 and 15 to 18), the Kantorovich averaging 
method (refs. 19 and 20), various perturbation methods (ref. 21 to 24), and 
incremental harmonic balance method (ref. 25). Studies based on the simpli- 
fied Berger's hypothesis (ref. 26) have also been made with the use of 
the Galerkin method (refs. 27 and 28). Yamaki et al. (ref. 18), Chia (ref. 
29) and Sythamoorthy (ref. 30) have presented their comprehensive reviews on 
both free and forced nonlinear vibrations of plates. 

The finite element method has proven to be an extremely powerful tool 
for complex structures. Applications of the finite element method to large 
amplitude free vibrations of beams and rectangular plates was first present- 
ed by Mei (refs. 31 and 32). The inplane tensile force induced by the 
transverse deflection alone was assumed to be constant for each individual 
beam or plate element. Nonlinear frequencies of beams and rectangular 

plates with various boundary conditions agreed well with the approximate 
continuum solutions of Woinowsky-Kr ieger (ref. 1), Chu and Herrmann (ref. 
3), Yamaki (ref. 15) and Evensen (ref. 33). Rao et al. proposed a novel 
scheme of linearizing the nonlinear strain-displacement relations in formu- 
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lating the nonlinear stiffness matrix. They studied nonlinear free vibra- 
tions of beams and circular plates (ref. 34) and rectangular plates (refs. 
35 and 36). Shear deformation and rotary inertia were also included in the 
formulation (refs. 37 and 38). Reddy and Stricklin (ref. 39) presented a 
linear and a quadratic isoparametric rectangular element using the linear- 
ized Reissner-type variational formulation to study large amplitude free 
plate vibrations. Inplane displacements were considered in their formu- 
lation. Two triangular elements have also been developed for nonlinear free 
vibrations of plates of arbitrary shape. The first one (ref. 40) is con- 
sistent with the higher-order bending element TRPLT1 (ref. 41) in NASTRAN®, 
and the second one (ref. 42) is consistent with the high-precision plate 
element of Cowper et al . (ref. 43). The solutions obtained for numerical 
examples include rectangular, circular, rhombic and isosceles triangular 
plates. Reddy and Chao (refs. 44 and 45) extended the earlier isoparametric 
rectangular elements to include transverse shear and laminated composite 

materials. Mei et al . (46) also extended the earlier triangular element to 

include laminated composite materials. 

Bhashyan and Prathap (ref. 47) and Sarma and Varadan (refs. 48 and 49) 
also presented a Galerkin and a Lagrange-type finite element formulation for 
nonlinear free vibrations of beams with ends restrained from longitudinal 
movement. They obtained frequency values at the instant of maximum ampli- 
tude which was based on a new criterion for defining nonlinear frequency 

presented in references 50 and 51. However, the results (refs. 47 to 50) do 
not agree with those classic continuum solutions (refs. 1 and 33). What 
they actually solved is a linear beam vibration problem subjected to an 
initial axial tensile force as commented on in reference 52. 

In this report, a finite element formulation is presented for large 
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amplitude vibrations of slender beams and thin plates subjected to harmonic 
excitation. Harmonic force matrices are developed for nonlinear vibrations 
of a beam and a rectangular plate element under uniform harmonic loading. 
Longitudinal or inplane deformation and inertia are both included in the 
formulation. These effects were neglected in the earlier finite element 
nonlinear free vibration analyses (refs. 31, 32, 34 to 38, 40, 42 and 46). 
Formulation of the harmonic force matrix [p] is based on the first order 
approximation of solutions of a Duffing system given by Hsu (ref. 53). He 
showed that the simple harmonic forcing function P q coswt is the first 
order approximation of the elliptic forcing function BAcn(A-r,n), and the 
well-known perturbation solution of a Duffing system to a simple harmonic 
forcing function is the first order approximation of the simple elliptic 
response Acn(AT,n). Derivation of the harmonic force and nonlinear stiff- 
ness matrices are given. Nonlinear forced response to uniform and simulated 
concentrated (uniform load on short beam segment) harmonic loadings, and 
improved nonlinear free vibration results are presented for beams and rec- 
tangular plates with various boundary conditions. Finite element results 
are compared with approximate solutions of elliptic function response, 
perturbation method and other numerical approaches. 
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SYMBOLS 


A 

B 


a, b 
a, b 

[C] 
c 

[D] 
d 

E 

(e> 

[F] 


fi, f 2 

[G] 

[H] 

h 

I 

[K] 

[K] 

[k] 

[k] 

L 

L 

i 

{M 0 } 


amplitude = w max /h for plate (w max /R for beam) 
nondimensional forcing amplitude factor 
length and width of rectangular plate 
length and width of rectangular plate element 
extensional material stiffness matrix 
constant defined in equations (46) and (47) 
bending material stiffness matrix 
length of the loaded element 
Young's modulus 
membrane strains 

linearizing function matrix defined in equation (27) 
linearizing functions 

matrix relating generalized coordinates and membrane 
strains defined in equation (26) 

matrix relating generalized coordinates and curvatures 
defined in equation (23) 

thickness of plate 

area moment of inertia of beam 

system linear stiffness matrix 

system nonlinear stiffness matrix 
element stiffness matrix 

element nonlinear stiffness matrix 
length of beam 

mathematical operator defined in equations (34) and (35) 
length of a beam element 
bending and twisting moments 
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[M] 

[m] 

{ N} 

[P] 

P o 

[p] 

[Q] 

R 

S 

[T b ], LTJ 

T 

t 

U 

u, v, w 
X, y, z 

a i* B i 
8 

Y 

{«> 

U) 

n 

{<} 

X 

V 

5 

P 

T 


system mass matrix 
element consistent mass matrix 
membrane resultant forces 
system harmonic force matrix 
nondimens ion al forcing paraneter 
element harmonic force matrix 

matrix relating generalized coordinates and linearizing 
function defined in equation (25) 

radius of gyration of beam cross-section 

beam cross-sectional area 

matrix relating element nodal displacements and generalized 
coordinates defined in equations (21) and (22) 

kinetic energy 

time 

strain energy 
displ acements 
coordinates 

general i zed coordi nates 

nonlinear coefficient in Duffing equation 

shear strain 

element nodal displacements 
strains 

modulus of elliptic function 
curvatures 

circular frequency of elliptic function 

Poisson's ratio 

norm 

mass density 
nondimensional time 
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4 > 

* 

0) 

Subscripts : 

b 

L 


mode shape normalized to unit value of the largest 
component 

stress function 

frequency 


bending 

linear 


m 


membrane 


FINITE ELEMENT FORMULATION 
Strain and Kinetic Energies 

From von Kantian' s large deflection plate theory, the nonlinear strain- 
displacement relations are defined as 


{e} = (el + z {<) 


( 1 ) 


given by 


midsurface strains 

(el 


"3_u 

+ I 

3w 2 
(— ) 


3 x 

2 

3X 

{e} =< 

>1 >> 
ro 1 <d 

H- I 
2 

3 W 2 
(— ) 
3y 


3_U 

+ 




3 X 

3 X 



r 3 2 W 




" 3X 2 


{<} = < 

3 2 W 

" 3/ 

!■ 


( 2 ) 


( 3 ) 


a 2 w 


3x3y J 


and u, v, w are displacements of the plate midsurface in the x, y, z - 
directions, respectively. 

The membrane or inplane resultant forces {N} and bending and twisting 
moments {M 0 } are related to the strains and curvatures by 
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= [C] (e} 


( 4 ) 


{N} = 


{M > 
o 


xy' 


My 


xy , 


CD] M 


( 5 ) 


where [C] and [D] are the extensional and bending material stiffness 
matrices, respectively. For plates of isotropic material and uniform thick- 
ness h, the material stiffness matrices are 


[C] - 


Eh 


l-v< 


1 v 0 
v 1 0 

1-v 


0 0 


2 J 


( 6 ) 


[D] - 


Eh 3 


12 (1-v 2 ) 


1 v 0 
v 1 0 

1-v 
2 


0 0 


( 7 ) 


where E and v are the Young's modulus and Poisson's ratio, respectively. 
The total strain energy for a plate element is given by 


u = u. + u 

b m 


( 8 ) 


with 


9 



= — // {M 0 ) T {<} dx dy 

= I ff {k} T [D] {<} dx dy (9) 

2 

U l = — ff {N}"*" { e> dx dy 
m 2 

= - ff (e> T [C] (e> dx dy (10) 

2 


where the subscripts b and m denote the bending and membrane components, 
respectively. 

The kinetic energy of a plate element executing harmonic oscillations 
is 


T = -L ph ff (u 2 + v 2 + w 2 ) dx dy (11) 

2 

where p is the mass density and (•) means differentiation with respect 
to time t. 

Similar expressions for a beam element can be obtained from the corre- 
sponding plate equations (1) to (11) by letting 

v = 0 


3y 
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[C] = ES 


[D] = El (12) 

where S and I are the area and moment of inertia of the beam cross- 
section, respectively. These beam equations are given in Appendix A. 

PI ate El ement 

The finite element used in the present formulation shown in figure 1 is 
the conforming rectangular plate with 24 degrees-of- freedom due to Bogner et 
al. (ref. 54). The displacement functions within the element are assumed as 

W = + a 2 X + a 3 y + a^x 2 + a 5 xy + a 8 y 2 

+ a 7 x 3 + a 8 x 2 y + agxy 2 + a 10 y 3 

+ anX 3 y + a^x 2 / + ai3Xy 3 

+ 014X3/ +a 15 x 2 y 3 + a 16 x 3 y 3 (13) 

u = + e 2 x + 3 3 y + Bi+xy (14) 

v = 6 5 + e 6 x + e 7 y + Bsxy (15) 


The twenty four generalized coordinates 




(16) 


{a} = [aj , a2 » • • • » ai6 ^ 

{0} T = [Pi, g 2 , .... 0 8 3 (17) 


can be determined from the element nodal displacements 



(«} T = [iy T , (<y T ] 

(18) 

with 

(5 b } T = [h 1# w 2> w 3 , w 4 , w xl , 



W yl’ W xyl’ W xy4^ 

(19) 


(V T = [Ui, Uz, U3. U*. v lf V 2 , V 3 , V 4 ] 

(20) 

The relationship 
displacements can 

between the generalized coordinates and 
be written as 

the element nodal 


{«) » [T b ] <«„} 

(21) 


( B ) - [T m ] «„> 

(22) 

Substituting 

equation (13) into equation (3), the 

curvatures can be 


expressed in terms of the generalized coordinates as 
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M = [H] {a} 


(23) 


The bending strain energy U b of equation (9) leads to the element linear 
bending stiffness matrix 


[k b ] - [T b ] T ff [H] T [D][H] dx dy [T b ] (24) 


Similarly, the kinetic energy T of equation (11) leads to the element 
consistent mass matrices [m b ] and [m m ], These three matrices can be 
found explicitly in reference 54. 

The membrane strain energy of equation (10) can be linearized 

utilizing the linearizing functions 



■ 1 M [T b ]{« b ! 


(25) 


where the element displacements { 5 b > are obtained from the plate deflec- 
1 tion through an iterative procedure discussed in the next section. Thus, 

I 

; the membrane strains of equation (2) become 

i 
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with 


r ^ 



9w 


9_U 



9x 


9x 



» CF]' 


1 +< 


9w 


sy 


lay 1 


9_U + 9v 

9 y ax 





• [F][q] («} * 

[G] 

- 1 [G]J 

1 

ko}\ 

m 


(26) 


[F] 


fl 0 

0 f 2 

f 2 f x 


(27) 


and the linearized membrane strain energy in terms of the element nodal 
displacements is 


4 [{5 » )T(5 " )Tl 


^ 0 


+ 


0 

0 


M 


<v| 

<y| 


( 28 ) 


where the element linear membrane stiffness matrix is given by 
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[k m ] - [T m ] T //[G] T [C][G]dx dy[T m ] (29) 


This [k^] matrix is also given explicitly in reference 54. The sub 
matrices of the linearized nonlinear stiffness matrix [k] are given by 

U b ] = [T b ] T // [q] T [F] T [C] [F] CQ] dx dy [T b ] 

(30) 

fftj * [T d ]T " W]T [F]T [C] [G] dx d * 

(31) 


(32) 


Evaluation of [k] is based on numerical integration using a four-point 
Gaussian integration which can exactly integrate for polynomial of cubic 
order. The matrices [T^], [T ], [H], [Q] and [G] defined in equations (21), 
(22), (23), (25) and (26) are given in Appendix B. 

The corresponding equations and matrices for the beam element are given 
in Appendix A. Element linear bending stiffness matrix [k^] , membrane 
stiffness matrix [k m ], consistent bending mass matrix [m^] ana membrane 
mass matrix [m m ] for the beam can be found in reference 55. 

Element Harmonic Force Matrix 

In classic continuum approach, the dynamic von Karman plate equations 
of motion are (refs. 2 and 3) 


^ = E (w, 2 xy 


w. 


XX 


W, 


yy 


) 


(33) 
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( 34 ) 


r(w,!j>) = phw, tt + DV'+W - h (*, w, xx 

+ ^’xx w *yy “ ^’xy w ’xy^ " ^ 

The equation of motion for a beam with ends restrained from longitudinal 
movement is (ref. 4) 

L(w) = P S w, u + El w, xxxx - Nw, xx - F(t) = 0 (35) 

with 


N = — f w,* dx (36) 

2L x 

For single mode approximate solutions, the deflection function is assuned in 
the form 

w = (hq(t) i»(x,y) , for plate (37a) 

1 Rq(t)<|)(x) , for beam (37 b) 

where R is the radius of gyration of beam cross-section, and mode shape <|> 
satisfies the related boundary conditions. Substitution of equation (37a) 
into equation (33), the stress function i|>(x,y) is then solved from the 
compatibility equation (33). Application of the Galerkin procedure 
//L(w,i|> )<|> (x,y)dxdy = 0 (or / L(w)<f>(x)dx = 0 for beam) yields a modal 
equation in the Duffing form (refs. 4, 5 and 15 to 17) 

mq, tt + kq + k q 3 = F(t) (38) 

or in nondimensional time x (ref. 53) 
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q» TT + q + eq 3 = F(t) 


(39) 


When the forcing function F(t) is simple harmonic P Q cos on, an approxi- 
mate solution of equation (39) using the perturbation method is the well- 
known result (refs. 4, 15, 53). 


(— ) 2 = 1+1 BA 2 - J 2 . (40) 
u> L 4 A 

With a simple elliptic forcing function F(t) = BAcn (Ax,n) = Bq as the 
external excitation to the Duffing system, an elliptic response (refs. 4, 5, 
15 to 17, and 53) 

q = A cn (At , n ) (41) 

is obtained as an exact solution of equation (39), where B is the non- 
dimensional forcing amplitude factor, A and n are the circular frequency 
and the modulus of the Jacobian elliptic function, and A = w m , /h (w m /R 
for beam) is the amplitude. By expanding the elliptic forcing function 
into the Fourier series and comparing the orders of magnitude of the various 
harmonic components, Hsu (ref. 53) showed that the single harmonic forcing 
function P q cos ojt is the first order approximation of the elliptic forcing 
function BAcn (At,n). Also the first order approximation of the elliptic 
response of equation (39) yields the same frequency-amplitude relations of 
equation (40) as the perturbation solution. In obtaining the exact elliptic 
response of equation (39), the excitation force F(x) = Bq is treated as a 
linear spring in the Duffing equation 
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(42) 


q, xT + (l-B)q + Bq 3 = 0 


This linear spring force Bq posses a potential energy of V = Bq 2 /2 . The 
potential energy of an element subjected to a uniform harmonic forcing 
function can thus be approximated by 


B 

— // w 2 dxdy , for pi ate 

2 

1 — / w 2 dx , for beam 

2 


(43a) 

(43b) 


An examination of equations (11) and (43) indicates that the harmonic force 
matrix for the plate or beam element under uniform loading F Q cosa)t is 


cF 


[p] = 


0 

Aph 2 

[m b ] 

, for pi ate 

(44a) 

cF 

0 

ARpS 

[m b ] 

, for beam 

(44b) 


The actual applied distributed force intensity F^N/m 2 or psi for plate, 
N/m or lb/in. for beam) is related to the dimensionless forcing parameter 
P q and the dimensionless forcing amplitude factor B by 


B 



A 


cF 

o_ 

Aph 2 ^ 2 


ARp S<u 2 


, for pi ate 


, for beam 


(45 a) 


(45b) 
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where 


c = — (L° a ded elements)^ dx( ty 
^ (Total plate area) 1 ** dxd ^ 

For plates under uniformly distributed force over entire plate 


c = // (Tota1 plate)* dxdy ( 47 j 

^ (Total plate area)** 1 dxdy 

which is simply the ratio of volumes (area) under plate (beam) mode shape 
and square of mode shape. The harmonic force matrix of equation (44) de- 
pends on the plate amplitude A = w max /h ( W max / R for beam) and P Q (or 

F o>* 

Equation of Motion and Solution Procedures 
The application of the Lagrange's equation leads to the equation of 
motion for the present rectangular plate and beam elements under the influ- 
ences of inertia, elastic, large deflection and uniform harmonic excitation 
force as 



0 

•• ( 

'[k b ] 

1 

o 



{6} + l 



0 

1 — 1 
3 
3 

V 

_ 0 



i — t 

£ k bJ 


1 1 

OL 
1 1 

i 

O 


0 


0 

L 

\ 

oj 


[ 6 } = 0 


(48) 





The coupling between bending and membrane stretching is evident by the pres- 
ence of [k bm J and [k mt) ] matrices in equation (48). Nonlinear free vibra- 
tion is a special case of the general nonlinear forced vibration problem 
with P Q or [p] = 0 in equation (48). 

By assembling the finite elements and applying the kinematic boundary 
conditions, the equations of motion for the linear free vibration of a given 
plate or beam may be written as 

“> L 2 CM] 4> 0 = [K] W Q (49) 

where [M] and [K] denote the system mass and linear stiffness matrices, 
respectively, u> L is the fundamental linear frequency, and denotes 
the corresponding linear mode shape normalized to unit value of the largest 
component. The deflection w max ^ 0 is then used to obtain the element 
nonlinear stiffness matrix [k] through equations (25 and 30 to 32). The 
element harmonic force matrix is obtained through equation (44) for given 
P Q and A. The nonlinear forced vibration response is approximated by a 
linearized eigenvalue equation of the form 

O) 2 [M] Wx = ([K] + [K] - [P]) {*}! (50) 

where a> is the fundamental nonlinear frequency associated with dimension- 
less amplitude A and force parameter P Q , and {<(>} i is the corresponding 
normalized mode shape of the first iteration. The iterative process can now 
be repeated with w max ^i until a convergence criterion is satisfied. 
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Convergence and Strains 


The three displacement convergence criteria proposed by Bergan and 
Clough (ref. 56) for nonlinear static and post-buckling analyses are employ- 
ed in the nonlinear free and forced vibrational problems. The three dis- 
placement norms are the modified absolute norm, modified Euclidean norm and 
the maximum norm. They are defined as 


and 


l« 


'abs 


i : 


A<t>. 


n 


1-1 


^i ,ref 


1 I 'euc 


"l n 

— z 

A<t> i 

1 

CNl 

n i=l 

^i ,ref 

J 


(51) 


(52) 


I Climax 


max 

i 


A<j> i 

^i ,ref 


(53) 


where n is the total degrees-of-freedom of the beam or plate structure, 
A<j> ^ is the change in displacement component i during interation cycle. 
The reference displacement <|>.j is the largest displacement component of 

the corresponding type (ref. 56). In addition a frequency norm is also 
introduced in the present study and it is defined as 



Ago 

U) 


(54) 
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where Aoj is the change in nonlinear frequency during the iteration cycle. 
A typical plot of the four norms versus nunber of iterations for a simply 
supported beam of slenderness ratio L/R = 100 with immovable end supports 
subjected to a uniform harmonic force of P Q = 2.0 at A = w max /R = 4.0 is 
shown in figure 2. Figure 3 shows the convergence characteristics of a 
simply supported square plate of 1 ength-to- thickness ratio a/h = 240 with 
immovable inplane edges (u = 0 at x = 0 and a, v = 0 at y = 0 and a) sub- 
jected to a uniform harmonic force of forcing parameter P Q = 0.2 at A = 

w /h = 1.0. All four norms exhibit the important behavior of straightness 
max 

and parallelism as described in reference 56. Therefore, there is no great 
significance as to what specific norm is being used since they are parallel, 
and also an upper bound or maximum error on displacement and frequency con- 
vergences can be estimated since they are straight. The results presented 
in this report, convergence is considered achieved whenever anyone of the 
norms reaches a value of 10 -5 . 

Once the convergence is satisfied, the membrane and bending strain 
components at the four nodes of the plate element can be obtained from the 
plate deflection and equations (2), (3), (21) to (23) and (25). They are 
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Maximum normal strain is evaluated with z at the top or bottom surfaces 
(z = ± h/2). The corresponding strain expressions for the beam element are 
given in Appendix A. 


RESULTS AND DISCUSSION 
Beams 

Improved Nonlinear Free Vibration 

The fundamental frequency ratios <o/u^ of free vibration at various 
amplitude A = W max /R 'for a simply supported beam (L/R = 100) with both 
ends immovable (u * 0 at x=0 and L) are shown in table I. Due to symmetry, 
only one-half of the beam which is divided equally into six elements is 
used. Finite element results with and without considering longitudinal 
deformation and inertia (LOI) in the analyses are both given. It shows that 
the effects of improved finite element results by including LDI in the 
formulation are to reduce the nonlinearity. The elliptic function solution 
(refs. 1, 4 and 53) is also given to demonstrate the closeness of the 
earlier finite element results without LDI. Raju et al. (ref. 57) used the 
Rayleigh-Ritz method in investigation of the effects of inplane deformation 
and inertia on large amplitude flexual vibration of slender beams. 
Appropriate frequency-amplitude relationship using Rayleigh-Ritz method is 
also given in table I. In Rayleigh-Ritz method, accuracy of the frequency 
depends heavily on the assumed beam deflection functions (usually linear 
mode shapes). The present finite element approach, the beam deflection 
approachs the "true" nonlinear beam deflection through an iterative process. 
The final beam deflection, therefore, would be more accurate, and also the 
frequency ratios obtained. 
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TABLE I. FREE VIBRATION FREQUENCY RATIOS «/« FOR A SIMPLY SUPPORTED BEAM 
WITH IMMOVABLE AXIAL ENDS. L 
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Nonlinear Forced Response of Beams with Immovable Ends 


Table II shows the frequency ratios of the same simply supported beam 
(L/R = 100) subjected to a uniform harmonic force of P Q = 2.0. It demon- 
strates the closeness between the earlier finite element results without 
LDI, the simple elliptic response (refs. 4 and 53) and the perturbation 
solution of equation (40). The present improved finite element results 
indicate clearly that the effects of LDI are to reduce the nonlinearity as 
in the free vibration case. Similar results of a clamped beam (L/R = 100) 
with immovable axial end supports subjected to a uniform harmonic force of 
P q =1.0 are also given in table II. The present finite element results of 
a slender beam (L/R = 100) to uniform harmonic excitation of P Q =0, 1.0 
and 2.0 are given in figures 4 and 5 for simply supported and clamped 
boundary conditions, respectively. 

Table III shows the maximum strain-ampl itude relation for the simply 
support and clamped beams (L/R = 100) with immovable axial end supports. 
These results are also shown in figure 6 for absolute value. The service 
life can thus be estimated from the maximum strain and frequency in conjunc- 
tion with the corresponding material strain-cycle fatigue data (S-N curve). 
This is not possible if small deflection free vibration analysis is employ- 
ed. 


Nonlinear Forced Response of Beams with Movable End-Support 

Figures 7(a) and (b) show the nondimensional amplitude A versus oi/u^ 
for a simply supported beam of slenderness ratio L/R = 100 and 20, respec- 
tively. One of the end supports (x = L) is assumed to be free to move in 
the axial direction. For a highly slender beam (L/R > 100), tne hard spring 
type nonlinearity due to large deflection is so small as shown in figure 
7(a), therefore, it can be practically neglected. Longitudinal deformation 
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Table II. FORCED VIBRATION FREQUENCY RATIOS w/u> FOR A SIMPLY SUPPORTED 
AND A CLAMPED BEAM WITH IMMOVABLE AXIAt ENDS. 



Without LDI a 

With LDI 


Simple 


Earlier 


Present 

w 

El liptic 


, Finite Element 

Finite Element 

. max 

Response 

Perturbation 

First 

Fi nal 

Final 

R 

(refs. 4 

Solution 



Result 




Iteration 

Result 



and 53) 





Simply Supporter 

— «- 

i Beam Subjected 

to P 0 = 2.0 (F ( 

5 = 2.31506 

x 10 5 N/m) 

- 1.0 

1.7852 

1.7854 

1.7852 

1.7856 

1 . 7682 ( 3 ) b 

± 2.0 

0.8472 

0.8660 

0.8621 

0.8460 

0.7108(4) 


1.6557 

1.6583 

1.6563 

1.6512 

1.5829(4) 

± 3.0 

1.4003 

1.4216 

1.4102 

1.3760 

1.2123(4) 


1.8217 

1.8314 

1.8226 

1.8002 

1.6743(4) 

± 4.0 

1.8413 

1.8708 

1.8453 

1.7846 

1.5871(6) 


2.1013 

2.1213 

2.0988 

2.0495 

1.8759(6) 

± 5.0 

2.2606 

2.2995 

2.2525 

2.1619 

1.9371(7) 


2.4361 

2.4673 

2.4236 

2.3432 

2.1337(7) 


Clamped Beam Subjected to P c 

, = 1.0 (F 0 = 5. 

,73862 x 10 5 N/m) 

± 1.0 

0.2118 

0.2165 

0.2096 

0.2091 

0.1772(3) 


1.4307 

1.4307 

1.4297 

1.4297 

1.4251(3) 

± 2.0 

0.8279 

0.8292 

0.8215 

0.8203 

0.7905(4) 


1.2987 

1.2990 

1.2942 

1.2936 

1.2743(4) 

± 3.0 

1.0401 

1.0433 

1.0279 

1.0239 

0.9726(5) 


1.3232 

1.3248 

1.3127 

1.3099 

1.2694(5) 

± 4.0 

1.2183 

1.2247 

1.1979 

1.1888 

1.1151(6) 


1.4101 

1.4142 

1.3910 

1.3836 

1.3197(6) 

± 5.0 

1.3938 

1.4042 

1.3619 

1.3457 

1.2513(8) 


1.5322 

1.5401 

1.5016 

1.4874 

1.4014(8) 


a. Longitudinal deformation and inertia. 

b. Number inside parenthesis denotes the number of iterations to get a converged 
solution. 


29 





















P 0 COSujr 


i 

I 

i 

I 


i 


I 

I 

I 


i 

! 

I 


OJ 



nx» aonindwv 


30 


Figure 4. Amplitude versus frequency for a simply supported beam with immovable axial end 
supports at uniform P = 0, 1.0 and 2.0. 



Table III. MAXIMUM STRAIN FOR BEAMS WITH IMMOVABLE AXIAL ENOS 

(L/R = 100) INCLUDING INPLANE DEFORMATION AND INERTIA. 


a = 

R 

Simply-Support 

Clamped 

1.0 

1.955 x 10‘3 

-4.741 x 10-3 

2.0 

4.386 x 10-3 

-9.367 x 10-3 

3.0 

7.292 x 10-3 

-14.107 x 10-3 

4.0 

10.676 x 10-3 

-19.169 x 10-3 

5.0 

14.535 x 10-3 

-24.726 x 10-3 
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and inertia effects are more pronounced in a short beam than a long one. As 
a consequence, the reduction of nonlinearity due to longitudinal deformation 
and inertia, from the nearly small deflection linear case in figure 7(a), 
leads to a situation that the beam eventually exhibits slightly soft spring 
type nonlinearity as shown in figure 7(b). Atluri (ref. 58) also obtained 
similar nonlinearity of softening type in his investigation. 

Concentrated Harmonic Force 

The element harmonic force matrix [p] derived in equation (44) is for 
uniformly distributed load over the element. The method used here to simu- 
late a concentrated force is to let the length of the loaded beam element 
become smaller and smaller. This is demonstrated by a concentrated force 
applied at the center of a simply supported beam (L/R = 100) with immovable 
inplane edges. The magnitude of the concentrated force is equal to the same 

beam under a uniform distributed load of P = 1.0 (F = 1.158 x 10 5 N/m or 

o o 

661 lb/in.) over the entire beam. Therefore, the uniform loading of the 

loaded element for the concentrated case is F = 115.8 x L/d kN/m where d 

o 

is the length of the loaded beam element. The constant c is given in 
equation (46). The results are given in table IV at various d/L ratios and 
in figure 8 for d/L = 5%. The elliptic function and perturbation solutions 
(without inplane deformation and inertia) are also given in table IV. It is 
shown that the concentrated force case is approximately 1.6 times as much 
severe as the uniform distributed force for the case studied. 
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Table IV. CONVERGENCE OF FREQUENCY RATIOS «/« WITH LOADED LENGTH d FOR A 

SIMPLY SUPPORTED BEAM (L/R = 100) WiVh IMMOVABLE INPLANE EDGES SUBJECTED 

TO A CONCENTRATED FORCE F Q = 115.8 x L/d kN/m AT THE CENTER. 


A = w max 

Withoui 

LDI* 

Finite Element with LDI 

Elliptic 

Solution 

Perturbation 

Solution 


at (c 

/L)% 


R 

20 

5 

1 

0.5 

- 1.0 

1.6607 

1.6608 

1.6316 

1.6411 

1.6423 

1.6425 

- 2.0 

.9695 

.9821 

.8585 

.8506 

.8497 

.8497 


1.5894 

1.5923 

1.5077 

1.5134 

1.5143 

1.5143 

- 3.0 

1.4519 

1.4710 

1.2742 

1.2718 

1.2717 

1.2717 


1.7815 

1.7920 

1.6273 

1.6317 

1.6325 

1.6326 

- 4.0 

1.8711 

1.8993 

1.6229 

1.6226 

1.6229 

1.6229 


2.0751 

2.0959 

1.8446 

1.8485 

1.8494 

1.8495 

5.0 

2.2801 

2.3181 

1.9606 

1.9615 

1.9620 

1.9621 


2.4179 

2.4498 

2.1116 

2.1154 

2.1164 

2.1165 


♦Longitudinal deformation and inertia. 
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Table VI. CONVERGENCE OF FREQUENCY RATIOS oj/oj WITH GRIDWORK REFINEMENT FOR 

A SIMPLY SUPPORT SQUARE PLATE ( a/h = 240) WITH IMMOVABLE 
INPLANE Edge Subjected to P Q = 0.2. 




Gridwork 


_ w max 
h 

2x2 

(4 Elements) 

3x3 

(9 Elements) 

4x4 

(16 Elements) 





± 0.2 

0.1645(3)* 

1.4248(3) 

0.1643(3) 

1.4238(3) 

0.1636(3) 

1.4237(3) 

± 0.4 

0.7815(3) 

1.2697(3) 

0.7800(3) 

1.2682(3) 

0.7792(3) 

1.2677(3) 

4 0.6 

0.9576(4) 

1.2588(4) 

0.9544(4) 

1.2560(4) 

0.9530(4) 

1.2550(4) 

± 0.8 

1.0937(5) 

1.3026(5) 

1.0886(5) 

1.2981(5) 

1.0865(5) 

1.2963(5) 

o 

• 

H 

+1 

1.2242(5) 

1.3781(5) 

1.2171(6) 

1.3717(6) 

1.2143(5) 

1.3691(5) 


*Number in parenthesis denotes the number of iterations to get a 
congerged solution. 
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Plates 


Improved Nonlinear Free Vibration 

The fundamental frequency ratios u>/a> L of free vibration at various 
amplitude A = w max /h for simply supported square (a/h = 240) and rectangu- 
lar (a/b = 2 and a/h = 480) plates with immovable inplane edges (u=0 at x = 
0 and a, v=0 at y=0 and b) are shown in table V. Due to symmetry, only one 
quarter of the plate modelled with 9 (or 3x3 gridwork) elements of equal 
sizes is used. Both finite element results with and without inplane defor- 
mation and inertia (IDI) are given. It shows that the improved finite ele- 
ment results by including IDI in the formulation are to reduce the non- 
linearity. The elliptic function solution and perturbation solution (with 
inplane deformation only, refs. 4 and 53) are also given to demonstrate the 
closeness of the earlier finite element results without IDI. Raju et al . 
(ref. 57) used the Rayl eigh-Ri tz method in investigation of the effects of 
IDI on large amplitude free flexural vibration of thin plates. The linear 
mode shape is very close to the nonlinear mode shape for the simply support- 
ed case. Therefore, the Rayl eigh-Ri tz solution demonstrates a good result 
compared to the present improved finite element solution. 

Convergence with Gridwork Refinement 

Table VI shows the frequency ratios for a simply supported square plate 
(a/h = 240) with immovable inplane edges subjected to a uniform harmonic 
force of P Q =0.2 with three finite element gridwork refinements. Only 
one quarter of the plate was used in the analysis due to synmetry. Exami- 
nation of the results shows that the present finite element formulation 
exhibits excellent convergence characteristics. Therefore, a 3 x 3 (or 9 
elements) in a quarter of plate was used in modeling the plates in the 
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Table V. FREE VIBRATION FREQUENCY RATIOS w/u> FOR A SIMPLY SUPPORTED PLATE 
WITH IMMOVABLE INPLANE EDGES. L 



Without IDI a 

With Inplane 
Deformation 
(No Inertia) 

With IDI 



Elliptic 




n _ w max 


Function 


Rayleigh 

Present 

h 

Finite 

Result 


Ritz 

Finite 


Element 

(refs. 4 

Perturbation 

Result 

Element 


Result 

and 53) 

Solution 

(ref. 57) 

Result 


Square Plate (a/h = 240) 

0.2 

1.0185 ( 3) b 

1 

1.0196 

1.0149 

1.0134(3) 

0.4 

1.0716(3) 


1.0761 

1.0583 

1.0528(3) 

0.6 

1.1533(4) 

1.1625 

1.1642 

1.1270 

1.1154(4) 

0.8 

1.2565(6) 

1.2734 

1.2774 

1.2166 

1.1979(5) 

1.0 

1.3752(7) 

1.4024 

1.4097 

1.3230 

1.2967(6) 


Rectangular Plate (a/b = 2, a/h 

= 480) 


0.2 

1.0238(3) 

1.0241 

1.0241 

1.0177 

1.0168(3) 

0.4 

1.0918(4) 

1.0927 

1.0933 

1.0690 

1.0658(4) 

0.6 

1.1957(6) 

1.1975 

1.1998 

1.1493 

1.1439(5) 

0.8 

1.3264(8) 

1.3293 

1.3347 

1.2533 

1.2467(6) 

1.0 

1.4762(11) 

1.4808 

1.4903 

1.3753 

1.3701(8) 


a. Inplane deformation and inertia. 

b. Number inside parenthesis denotes the number of iterations to get a 
converged solution. 
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remainder of the nonlinear forced responses presented unless otherwise 
specified. 

Nonlinear Forced Response of Plates with Immovable Inplane Edges 

Table VII shows the frequency ratios <o/u>^ for simply supported and 

clamped square plates (a/h = 240) subjected to a uniform harmonic force of 

P =0.2. It demonstrates the closeness between the earlier finite element 
o 

formulation without IDI, the simple elliptic response (refs. 4 and 53) and 
the perturbation solution (with inplane deformation only). The present 
improved finite element results indicate clearly that the effects of IDI are 
to reduce the nonl inearity. The present finite element results of a square 
plate (a/h = 240) to uniform harmonic excitation of P Q =0, 0.1 and 0.2 are 

given in figures 9 and 10 for simply supported and clamped boundary con- 

ditions, respectively. 

Nonlinear Forced Response of Plates with Movable Inplane Edges 

The dimensionless amplitude A versus the fundamental frequency ratio 
for a simply supported square plate (a/h = 240) with movable inplane 

edges subjected to uniform harmonic load P Q = 0, 0.1 and 0.2 is shown in 

figure 11. The nonlinearity is greatly reduced with the inplane edges no 
longer restrained as compared to the case of immovable inplane edges in 
figure 9. 

Concentrated Harmonic Force 

Application of the present finite element to the case of a concentrated 
force is to let the area of the loaded element becoming smaller and smaller. 
It is demonstrated by a concentrated force applied at the center of a simply 
supported square plate (a/h = 240) with immovable inplane edges. The magni- 
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Table VII. FORCED VIBRATION FREQUENCY RATIOS w/uj FOR A SQUARE PLATE 

(a/h = 240) WITH IMMOVABLE INPLANE EDGE^ SUBJECTED TO P Q = 0.2. 





Finite 

El ement 


Simple 





El liptic 

Perturbation 





Solution 

Without 

With 

_ w max 
h 

Response 
(refs. 4 


IDI a 

IDI 


and 53) 







Simply Supports 

id 


± 0.2 


0.1987 

0.1932(3) b 

0.1643(3) 


■ 

1.4281 

1.4274(3) 

1.4238(3) 

• 

o 

+1 


1 

0.8052(3) 

0.7800(3) 


1.2874 


1.2839(3) 

1.2682(3) 

± 0.6 

1.0084 

1.0110 

0.9984(4) 

0.9544(4) 


1.2983 

1.2995 

1.2898(4) 

1.2560(4) 

1+ 

o 

• 

00 

1.1703 

1.1755 

1.1528(6) 

1.0886(5) 


1.3686 

1.3718 

1.3524(6) 

1.2981(5) 

1+ 

• 

o 

1.3283 

1.3369 

1.3004(7) 

1.2171(6) 


1.4726 

1.4789 

1.4460(7) | 

1.3717(6) 



Cl amped 



± 0.2 

0.1200 


0.1180(2) 

0.1033(3) 


1.4195 

1.4195 

1.4195(2) 

1.4183(3) 

<3* 

• 

o 

+1 

0.7483 

0.7438 

0.7459(3) 

0.7372(4) 


1.2490 

1.2491 

1.2477(3) 

1.2426(4) 

1+ 

o 

• 

CTY 

0.8951 

0.8956 

0.8905(4) 

0.8746(4) 


1.2117 

1.2119 

1.2083(4) 

1.1966(4) 

± 0.8 

0.9941 

0.9954 

0.9863(5) 

0.9617(5) 


1.2203 

1.2210 

1.2137(5) 

1.1938(5) 

± 1.0 

1.0822 

1.0845 

1.0700(6) 

1.0362(5) 


1.2540 

1.2555 

1.2429(6) 

1.2140(5) 


a. Inplane deformation and inertia. 

b. Number inside parenthesis denotes the number of iterations to get a 
converged solution. 
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Figure 10. Amplitude versus frequency for a clamped square plate (a/h = 240) with immovable 
inplane edges at P n = 0, 0.1 and 0.2. 
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Figure 11. Amplitude versus frequency for a simply supported square plate (a/h = 240) with movable 
inplane edges at P = 0, 0.1 and 0.2. 


tude of the concentrated force is equal to the same plate under a uniformly 
distributed harmonic loading of = 0.1 (F = 45.74 N/m 2 or 0.66347X10- 2 
psi) over the total plate area. Therefore, the uniform loading of the load- 
ed element for the concentrated case is F q = 45. 74 (a/d) 2 N/n? where d 
is the length of the loaded square element. Table VIII gives the fundamen- 
tal frequency ratios co/uj l at (d/a) 2 = 16.0, 4.0, 1.0 and 0.25%. It 
indicates that the convergence is rapid and (d/a) 2 = 1.0% would yield accu- 
rate frequency response. Results obtained using earlier finite element 
without IDI and elliptic function (with in plane deformation but no inplane 
inertia) are also given. Nonlinear response of concentrated force obtained 
with (d/a) 2 = 1.0% is plotted in figure 12. Frequency ratios of the same 
plate to uniform harmonic force P Q =0.2 is also given. It shows that the 
concentrated force is approximately two to three times as much severe as the 
uniformly distributed force for the case studied. 

CONCLUSIONS 

The finite element method has been extended to analyze nonlinear forced 
vibration problems. Harmonic force matrices were developed for a beam and a 
rectangul ar plate element subjected to uniform harmonic excitation. 
Improved finite element results on nonlinear free flexural vibration of 
slender beams and thin plates are achieved by considering inplane deforma- 
tion and inertia effects in the formulation. Nonlinear free vibration can 
be simply treated as a limiting case of the more general forced vibration 
problan by setting the harmonic force matrix equal to zero. The effect of 
midplane stretching due to large deflection is to increase the nonlinearity, 
however, the effects of inplane deformation and inertia are to reduce non- 
1 inearity. 
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For beams with immovable end supports, only hardening type nonlinearity 
is observed. For beams of large slenderness ratio with a movable axial end 
support, the increase in nonlinearity due to large deflection is partially 
compensated by the reduction in nonlinearity due to longitudinal deformation 
and inertia. This leads to a negligible hardening type nonlinearity, and 
thus small deflection linear solution can be used, i.e. a simply supported 
beam with movable axial end support. For beams of small slenderness ratio, 
however, softening type nonlinearity is observed. 

For rectangular plates with immovable or movable inplane edges, how- 
ever, only hardening type nonlinearity is observed. This is because the 
middle surface of the bent plate is not a developable surface. 

For concentrated loading (beams and rectangular plates) yields re- 
sponses several times as severe as the uniformly distributed load. Fatigue 
life can be estimated with material fatigue S-N data. 
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APPENDIX A: BEAM ELEMENT 


The nonlinear strain-displacement relation is 


e = e + z < 


3u 1 /-SWso 3^ 

— + - (— ) 2 - Z — 

3 X 2 3 X 3X 2 


Axial resultant force and bending moment are given by 


(Al) 


N = ESe 


. ES [il + I 

3 X 2 3x 


(A2) 


M = El k 


= -El 


3 2 W 

3x 2 


(A3) 


The bending strain, membrane strain and kinetic energies are 


a 

U. = — / (— ) 2 dx (A4) 

b 2 0 3x 2 

U m / [l± +I(!l)2]2 dx (A5) 

2 0 3x 2 3 x 
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I 

t The element nodal displacements at the two end nodes are 

i 

i 

m T - [{«/ {y T ] 

I 

I 

= [w r W xl , W v W x2 , U v U 2 ] ( A9) 

The matrices relate the element nodal displacements and the generalized 
coordinates are 

10 0 0 

0 1 0 0 

_3 _ 2 _ l 

l 2 l l z l 

2 _ 1 _ _ _2 1 _ 

£ 3 £ 2 Z 3 Z 2 
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l o 

1 j, 

i % 


(All) 


The linearizing function is evaluated from the expression 


f a j. 9w 


2 9x 


- CQ]{«> 

2 


= - [0 1 2x 3x 2 ] [T b ]{6 } (A12) 


The matrix relates generalized coordinates and membrane strain defined in 
equation (26) is 


[a] ■= [o 

1] 

force matrix is 


' 156 


cF q * 22* 

4* 2 

42 OAR 54 

13* 

.-13* 

-3* 2 

the two end nodes 

of a 


(A13) 


symmetric 

156 

-22* 4* 2 


(A14) 


0 f -2z 0 0 

0 f 2(*f-z) 3(* 2 f-2* z) 0 


{6} (A15) 
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APPENDIX B: PLATE ELEMENT 


The inverse of matrix [T b ] in equation (21) is given by 


°1 







*8 




®1 2 




a i6 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 

a 

0 

p 

0 

0 

73 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 
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b 

7 2 

ab 

b 2 

I 3 

7 2 b 

ab 2 

b 3 

Pb 

7 2 P 

ab 3 

7 3 b 2 

7 2 b 3 

a 3 b 3 


1 

0 

b 

0 

0 

b 2 

0 

0 

0 

b 3 

0 

0 

0 

0 

0 

0 


0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

1 

0 

2a 

0 

0 

3a 2 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

1 

0 

2a 

b 

0 

3a 2 

2ab 

b 2 

0 

3a 2 b 

27b 2 

b 3 

3l 2 b 2 

27b 3 

37 2 b 3 


0 

1 

0 

0 

b 

0 

0 

0 

b 2 

0 

0 

0 
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0 

0 

0 
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0 

0 
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0 

0 

0 

0 
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0 

0 
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1 

0 

0 

2a 

0 

0 

3P 

0 

0 

0 
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1 

0 

0 
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(Bl) 

where 

a 

and 

F 

are 

the length 

and 

width 

Of 

rectangul ar pi ate 

el ement. 
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Matrix [H] in equation (23) is of the form 


| 



0 2 0 0 6x 

0 0 0 2 0 

0 0 2 0 0 


a 8 




°12 

2y 

0 

0 

6xy 

2X 2 

0 

2x 

6y 

0 

2x2 

4x 

4y 

0 

6x 2 

8xy 





a l6 

0 

6xy 2 

2y 3 

6xy 3 

6xy 

2x 3 

6xy 2 

6x 3 y 

6y 2 

12x 2 y 

12xy 2 

18x 2 y 2 


(B3) 


Matrix [Q] in equation (25) is 

a 8 a l 2 a l 6 

y 0 3x 2 2xy y 2 0 3x 2 y 2xy2 y 3 3x 2 y 2 2xy 3 3x 2 y 3 

x 2y 0 x 2 2xy 3y 2 x 3 2x 2 y 3xy 2 2x 3 y 3x 2 y 2 3x 3 y 2 

( B4) 
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Matrix [G] in equation (26) is given by 


Bl 64 38 

0 1 0 y 0 0 0 0 

0000001X 
OOlXOlOy 


i 





bb 
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